﻿namespace Bioinformatics

open DiffLCS

open System

/// Linear space global alignment
type LinearGlobalAlignment(a : DNABase[], b : DNABase[], subcost : int[,], gapcost : int) =
    // TODO: Share code with GlobalAlignment
    do
        if subcost.GetLength(0) <> 4 || subcost.GetLength(1) <> 4 then
            raise (new ArgumentException("The subcost matrix is not suitable."))

    let rec local a b =
        // Take the shortest as b to minimize memory consumption
        let a, b = if Array.length a < Array.length b
                   then b, a
                   else a, b
        let n, m = Array.length a, Array.length b

        if n = 1 || m = 1 then
            let nw = new GlobalAlignment(a, b, subcost, gapcost)
            nw.OptimalTrace()
        else
            let F = Array2D.zeroCreate 2 (m + 1)
            let c = Array2D.zeroCreate 2 (m + 1)
        
            let u = n / 2

            let mutable col' = 0
            let mutable col = 1

            // Writing to the row which will become the previous
            // in the first iteration in the double-loop.
            for j = 1 to m do
                F.[col, j] <- gapcost * j

            for i = 1 to n do
                // Exchange previous and current column
                col' <- col' ^^^ 1
                col <- col ^^^ 1

                F.[col, 0] <- gapcost * i

                for j = 1 to m do
                    // Match/substitution
                    let q = F.[col', j - 1] + subcost.[int a.[i - 1], int b.[j - 1]]
                    let p = F.[col', j] + gapcost // Deletion
                    let r = F.[col, j - 1] + gapcost // Insertion

                    let max0 = max q (max p r)
                    let i', j' = if max0 = q then col', j - 1
                                 elif p = max0 then col', j
                                 else col, j - 1

                    F.[col, j] <- max0

                    if i = u then
                        c.[col, j] <- j
                    else
                        c.[col, j] <- c.[i', j']

            let v = c.[col, m]

            (local a.[0 .. u - 1] b.[0 .. v - 1]) @ (local a.[u .. n - 1] b.[v .. m - 1])
            

    /// Returns the trace of an optimal alignment.
    member public this.OptimalTrace() =
        local a b

    member public this.Optimal () =
        let optimalalignment = this.OptimalTrace()
       
        let top = List.map fst optimalalignment
        let bottom = List.map snd optimalalignment

        let mapping = function
                        | -1 -> '-'
                        | x -> enum<DNABase>(x).ToString().[0]
        let f xs = xs |> List.map mapping |> List.toString

        (f top, f bottom)